A Circumplanetary Disk around PDS70c

Recent surveys have revealed that, almost ubiquitously, protoplanetary disks appear highly structured with rings and gaps, spiral arms, and asymmetries (e.g., Garufi et al. 2018; Andrews 2020). While other scenarios are discussed, these features are often interpreted as resulting from the presence of planets embedded in disks (e.g., Dong et al. 2015; Bae et al. 2018; Lodato et al. 2019). Additional observational support for such a scenario can be found in the form of local perturbation of the gas velocity field from Keplerian rotation (Pinte et al. 2018; Casassus & Pérez 2019; Teague et al. 2019). The quest to detect protoplanets embedded in their host disk through direct imaging has been challenging, with current detection limits on the order of a few Jupiter masses (M_Jup) at large radii (e.g., Huélamo et al. 2018; Asensio-Torres et al. 2021). A few protoplanet candidates have been claimed in the infrared (IR) and in the H_α line (e.g., Sallum et al. 2015; Reggiani et al. 2018) but remain controversial (Mendigutía et al. 2018).

The first robust detection through direct imaging techniques of a protoplanet still embedded in its natal disk was obtained in the young system PDS 70 (spectral type K7; M ∼ 0.8 M_⊙; age ∼ 5.4 Myr old; Müller et al. 2018) located at ∼112.4 pc (Gaia Collaboration et al. 2021) in the Upper Centaurus Lupus association (Pecaut & Mamajek 2016). PDS 70 b was discovered with an orbital radius of ∼22 au, and imaged at multiple IR wavelengths (Keppler et al. 2018; Müller et al. 2018) as well as in a filter centered on the Hα line (Wagner et al. 2018a). Subsequently, PDS 70 c was discovered in Hα imaging at the outer edge of the cavity with an orbital radius of ∼34 au (Haffert et al. 2019). These two planets carve a large cavity in the disk, evidenced by a cavity in dust (e.g., Dong et al. 2012; Hashimoto et al. 2012) and a gap in the ¹²CO gas emission along the orbit of PDS 70 b (Keppler et al. 2019) that indicates significant gas depletion. Observations and hydrodynamic simulations indicate that the planets' orbital configuration is stable, close to a 2:1 mean motion resonance, with PDS 70 b in a slightly eccentric orbit (e ∼ 0.2; Bae et al. 2019; Toci et al. 2020; Wang et al. 2021). The masses of the two planets are still uncertain, although both planets are likely lighter than 10 M_Jup to ensure dynamical stability (Wang et al. 2021) and a non-eccentric outer disk (Bae et al. 2019). Spectro-photometric analyses, limited to the IR regime (1–5 μm), remain inconclusive, but suggest planet masses between 1 and a few M_Jup (e.g., Müller et al. 2018; Mesa et al. 2019; Stolker et al. 2020) as well as a clear contribution from dust grains in clouds and/or circumplanetary disks (CPDs; Christiaens et al. 2019; Stolker et al. 2020; Wang et al. 2020).

CPDs play a fundamental role in planet formation, as they regulate the gas accretion onto the planet and determine the conditions for satellite formation. As gas enters the planet's sphere of influence, it falls at supersonic velocities onto the surface of the CPD (Tanigawa et al. 2012; Szulágyi & Mordasini 2017), possibly episodically (Gressel et al. 2013), leading to shocks that can ionize hydrogen and be traced in the Hα line. From observations of the Hα line, PDS 70 b and PDS 70 c are found to be accreting material from their host disk at a rate of ∼10⁻⁸ M_Jup per year (Wagner et al. 2018b; Aoyama & Ikoma 2019; Haffert et al. 2019; Thanathibodee et al. 2019; Hashimoto et al. 2020). Using Atacama Large Millimeter/submillimeter Array (ALMA) observations at ∼67 mas × 50 mas resolution, Isella et al. (2019) showed evidence for submillimeter continuum emission co-located with PDS 70 c, interpreted as tracing a dusty CPD, and for another compact continuum emission source located at ∼74 mas offset in a southwest direction from b. The emission around c, however, was not spatially separated from the outer ring. In this Letter, we present new ALMA observations with 20 mas resolution that provide an independent detection of a compact source of emission co-located with PDS 70 c and of low surface brightness emission within the cavity close to PDS 70 b. The Letter is organized as follows. Section 2 presents the observations and the procedure to calibrate the data. Section 3 presents our new images and analysis. Finally, we discuss our findings in Section 4.

This Letter presents new ALMA observations, hereafter referred to as LB19 ("Long Baselines 2019"), obtained in Band 7 (λ = 855 μm), under a Director's Discretionary Time program with ID 2018.A.00030.S. PDS 70 was observed during four execution blocks (EBs) with the C-8 configuration on 2019 July 27, 28, and 30, for a total on-source time of 43 minutes per execution. An observing log including the precipitable water vapor (PWV) levels and calibrator names is given in Appendix A.1. The spectral setup was tuned to optimize continuum detection, but includes the ¹²CO J = 3 − 2 line at 345.8 GHz and the HCO⁺ J = 4 − 3 line at 356.7 GHz, which will be presented in forthcoming papers. The raw data calibration was done with the CASA v.5.6.1 pipeline (McMullin et al. 2007) and the self-calibration and post-processing imaging were done using CASA v.5.4.0. We first flagged the channels that included the ¹²CO and the HCO⁺ lines and spectrally averaged the remaining channels to produce a continuum data set. We imaged the resulting visibilities with the tclean task using the multiscale CLEAN algorithm with scales of 0, 1, 3, and 6 times the beam FWHM, and an elliptic CLEAN mask encompassing the disk emission. To reduce the size of the data, we time averaged it to 6.06 seconds, i.e., 3 times the original integration time. After imaging, one EB image appeared of much lower signal-to-noise ratio (S/N), and therefore the corresponding visibilities were rejected. The individual images of the three remaining EBs (0, 1, and 3) did not appear astrometrically offset with respect to each other, which is as expected because they were taken very close in time. As the fluxes of all EBs match within 2%, we concatenated the three EBs and self-calibrated them all together. To determine a good initial model for the self-calibration, we used multiscale cleaning with the tclean task using a threshold of ∼7 times the rms noise level of the image. Using the tasks gaincal and applycal, we corrected for phase offsets between spectral windows, and between polarizations considering a solution interval of the scan length (solint = inf). Another iteration of phase self-calibration was done with a solution interval of 30 s. We reached an overall improvement in peak S/N of 34% after self-calibrating the LB19 data.

The LB19 data were combined with archival observations previously published in Isella et al. (2019) and are summarized in Table 1. These observations correspond to program ID 2015.1.00888.S (PI: E. Akiyama), taken in 2016 August and labeled SB16 ("Short Baselines 2016"), and to program ID 2017.A.00006.S (PI: M. Keppler) taken in 2017 December, labeled IB17 ("Intermediate Baselines 2017"). We refer the reader to Appendix A of Isella et al. (2019), where the procedure for the self-calibration of SB16 and IB17 data is described in detail. For all data sets, we used the statwt task to weight the visibilities according to their scatter. Before combining the LB19 data with the previously published data, we fitted an elliptical ring to the maximum of the outer ring in the image plane, for all data sets separately, to derive the center of the image. We then used the fixvis task to shift the image to the phase center and assign it to a common phase center using the fixplanets task on the center coordinate derived by Isella et al. (2019). The fluxes of the executions in LB19 differed by ∼3% from the archival data sets (IB17+SB16; Isella et al. 2019) and were rescaled using the DSHARP rescale_flux function. ¹⁹ After concatenation of the data, we followed the same procedure as explained above, with three rounds of phase self-calibration.

We proceeded with imaging of the final data using CLEAN. In a normal CLEAN workflow, after the CLEAN iterations terminate when the peak value of the residual image drops below a threshold value (4 × rms noise level in the observations considered here), a restored CLEAN model is combined with the residual image to form the CLEANed image. As discussed in Czekala et al. 2021; however, the units of these two quantities differ: the units of the restored CLEAN model are Jy {CLEAN beam}⁻¹ while the units of the residual image are Jy {dirty beam}⁻¹, as it originated as the dirty image. When the CLEAN beam (typically chosen to be an elliptical Gaussian) poorly approximates the dirty beam (as is common with multi-configuration ALMA data sets), the normal CLEAN workflow produces a CLEANed image with an incorrect flux scale and compromised image fidelity, especially for faint emission. This phenomenon was first described in Jorsater & van Moorsel (1995), and so we term it the "JvM effect". To correct for the unit mismatch, before combining the residual image with the restored CLEAN model, we first rescaled the residual image by the ratio of the CLEAN beam/dirty beam "volumes" (see "JvM correction", Czekala et al. 2021).

To test the effect of the angular resolution on the image features and assess their robustness, we performed a grid of CLEANed and JvM-corrected images, using Briggs weighting (Briggs & Cornwell 1992) with different robust parameters. A gallery of continuum images (and corresponding fluxes), synthesized from the new data set alone (LB19) and from data set combinations including the observations published by Isella et al. (2019) (IB17+SB16; LB19+IB17+SB16) is given in Appendix A.2. Depending on the data set and the robust parameter, our JvM-corrected images have a rms ranging between ∼4 and ∼26 μJy beam⁻¹ across beam sizes of 93 mas × 74 mas to 20 mas × 20 mas (Table 4 in Appendix A). We note that while the uv plane coverage and sensitivity are maximized when all data sets are combined (LB19+IB17+SB16), such a combination does not take into account the intrinsic changes of the emission that are due to the rotation of the system, and the change in the location of the dust surrounding the planets. Based on the orbital solutions of Wang et al. (2021), we expect a motion of ∼14 mas for both planets between 2017 December and 2019 July.

3.1. Continuum Images

Figure 1 presents a selection of images of the continuum emission of PDS 70 at 855 μm, synthesized from the new ALMA observations combined with short baseline data (LB19+SB16). The disk is well detected with a spatially integrated flux density of ∼176 ± 18 mJy (all images give similar values). After deprojecting the image with an inclination of ∼517 and a position angle of ∼1604 (Keppler et al. 2019), we computed an azimuthally averaged radial profile and found that the outer disk resolves in a ring extending radially from ∼04 (45 au) and ∼09 (100 au). The outer disk is not radially symmetric and shows a clear azimuthal asymmetric feature in the northwest (∼27% brighter at peak compared to the mean ring value), as already discussed by Long et al. (2018) and Keppler et al. (2019). When imaged at high resolution, the outer disk resolves into a narrow and bright ring with a faint inner shoulder detected in the image at the 3–4σ level (Appendix A.2). To better assess the presence of such substructures, we model the azimuthally averaged radial visibility profile using the frank package (Jennings et al. 2020). Our analysis, presented in Appendix B, recovers a double peaked profile for the outer disk. Such a substructure was already hinted in the data presented in Keppler et al. (2019). Inward of the outer disk, the dust-depleted cavity includes an inner disk that radially extends up to 016 (18 au) and presents faint additional emission in the west and in the south of the inner disk that will be discussed in the next subsection.

Zoom In Zoom Out Reset image size

3.2. Emission within the Cavity

Within the cavity, the inner disk appears well resolved with an integrated flux ranging between 727 ± 27 and 888 ± 59 μJy depending on the data set (Table 5 in Appendix A). When imaged at high angular resolution (e.g., Figure 1, left panel), it appears irregular and the emission is discontinuous in the north.

Continuum emission is also detected near the locations of the planets, confirming the findings of Isella et al. (2019). We use the same nomenclature as Isella et al. (2019) and label the continuum emission located close to planet b and c, b_smm and c_smm, respectively. The continuum emission around PDS 70 c, c_smm, is recovered in all images, and in particular in the new stand-alone high-resolution data set (LB19), where it appears as a 5.4 to 16σ feature depending on the robust parameter. c_smm clearly separates from the outer disk when imaged at resolutions finer than ∼40 mas. It appears unresolved even at our best angular resolution (∼20 mas; ∼2.3 au). We find that its peak intensity is similar in all the images that spatially resolve it from the outer disk (see Appendix A.3), confirming its point-source nature. Depending on the data set (IB17+SB16 or LB19+SB16) and the robust parameter, its peak intensity ranges between 80 ± 6 and 107 ± 15 μJybeam⁻¹. In the following, we will consider 86 ± 16 μJybeam⁻¹ as a reference for further discussion.

On the other hand, the emission located near PDS 70 b, b_smm, is only recovered when the new high-resolution data is combined with short baselines, and when the beam is larger than ∼50 mas. This indicates that it is low surface brightness, extended emission. Its peak intensity and morphology vary greatly between images of different data sets (Table 5 in Appendix A), which makes its morphology and properties difficult to recover accurately.

In order to assess whether the signal within the cavity could result from imaging artifacts, following Andrews et al. (2018) we subtracted the Fourier transform of the CLEAN model of the outer disk, after blanking out the pixels within the cavity (using an elliptical mask of 025 × 04), and image and model the visibilities carrying the residual signal from within the cavity. Figure 2 show two residual images, hereafter called "cavity images", for LB19+SB16 and IB17+SB16, that clearly show that the inner disk emission and c_smm are recovered in both epochs, the latter with a significance up to 18σ. On the other hand, b_smm is detected at a 3σ level only in some cavity images obtained from combined data sets. A gallery of cavity images is given in Appendix A.2.

Zoom In Zoom Out Reset image size

As an additional test, we perform a model fit of the cavity visibilities using the data set LB19+SB16, obtained after subtracting the Fourier transform of the CLEAN model of the outer disk using a robust parameter of 1. We consider a simple model for all three sources of emission within the cavity (namely the inner disk, b_smm, and c_smm), compute the Fourier transform using galario (Tazzari et al. 2018), and explore the parameter space using the Monte Carlo Markov chains implementation in emcee (Foreman-Mackey et al. 2013). Our model consists in a Gaussian ring for the inner disk, which enables the modeling of an additional structure within the inner disk, a point source for c_smm (between PA = 250° and 280°), and a circular Gaussian for b_smm, located in the south (between PA = 70° and 250°). A uniform prior was used over the allowed range for each parameter. Our best-fit model and residual maps are shown in Figure 3, and corresponding parameters are in Table 2. We find that the best-fit location of c_smm is Δ(R.A., decl. ) = ($-{215.1}_{-1.6}^{+1.8}$, ${37.8}_{-3.7}^{+3.3}$) mas, close to the predicted position of PDS 70 c Δ(R.A., decl. ) = (−214.8, 31.9) mas (see Appendix C). For b_smm, the location is constrained to Δ(R.A., decl. ) = (${70.1}_{-2.5}^{+2.4}$, $-{163.0}_{-3.2}^{+3.4}$) mas, offset from the predicted position of PDS 70 b (Δ(R.A., decl. ) = (96.9, −153.7) mas).

Zoom In Zoom Out Reset image size

Table 2. Best-fit Parameters for the Model to the Cavity Data for the Data Sets LB19+SB16 and IB17+SB16, with the 1σ Error

Data Set	finn	rinn	σinn	rb	θb	fb	rc	θc	fc
	(mJy)	(mas)	(mas)	(mas)	(deg)	(μJy)	(mas)	(deg)	(μJy)
LB19+SB16	${0.846}_{-0.047}^{+0.036}$	${2.0}_{-1.6}^{+25.0}$	${59.3}_{-12.0}^{+2.7}$	${178.5}_{-3.8}^{+2.7}$	${174.0}_{-1.3}^{+1.4}$	${83.1}_{-15.8}^{+12.4}$	${324.9}_{-2.7}^{+2.7}$	$-{70.6}_{-0.7}^{+0.6}$	${111.5}_{-13.6}^{+14.0}$
IB17+SB16	${0.765}_{-0.040}^{+0.018}$	${2.8}_{-1.0}^{+24.8}$	${52.6}_{-12.0}^{+0.1}$	⋯	⋯	⋯	${329.4}_{-10.1}^{+10.8}$	$-{68.9}_{-1.0}^{+1.1}$	${91.6}_{-13.1}^{+14.4}$
				Δ R.A. (mas)	Δ Decl. (mas)		Δ R.A. (mas)	Δ Decl. (mas)
LB19+SB16				${70.1}_{-2.5}^{+2.4}$	$-{163.0}_{-3.2}^{+3.4}$		$-{215.1}_{-1.6}^{+1.8}$	${37.8}_{-3.7}^{+3.3}$
IB17+SB16				⋯	⋯		$-{219.2}_{-6.5}^{+7.0}$	${47.9}_{-4.8}^{+4.9}$

Note. The flux, radial peak position, and width of the Gaussian for the inner disk are f_inn, r_inn, and σ_inn, respectively. The total flux and polar coordinates in the disk plane of b_smm and c_smm are f_b , r_b , and θ_b and f_c , σ_c , r_c , and θ_c , respectively. The relative apparent astrometry Δ(R. A., decl. ) is also provided.

Download table as:  ASCIITypeset image

From the orbital fits of Wang et al. (2021), the expected motions of the planets between the epoch of the long baselines observations (2017 December and 2019 July) is similar for both, ∼14 mas, which is smaller than the angular resolution of our observations. To search for possible motion of c_smm between the two epochs, we performed the same modeling as above on the IB17+SB16 data set. b_smm was not recovered in this fit, but the inner disk and c_smm were recovered. Using the best-fit positions for c_smm at the two epochs, and considering a 2 mas error in the centering of the two data sets, we find marginal evidence for a movement of the peak position of 10.9 ± 6.9 mas. We note that the nominal positional accuracy is defined as beam_FWHM/S/N/0.9 (Thompson et al. 2017, and ALMA Cycle 8 2021 Technical Handbook), with 0.9 a factor to account for a nominal 10% signal decorrelation. We consider two images in which c_smm is imaged at a decent S/N and separated from the outer disk, LB19+SB16 (r = 0.5) and IB17+SB16 (r = −0.3). With corresponding S/N of 8.9σ and 7.1σ on the peak intensity of c_smm, respectively, and a beam FWHM of 36 and 60 mas, respectively, the positional accuracies are ∼4.5 mas and 9.4 mas, respectively, comparable to the uncertainty that we derived for the apparent displacement of c_smm. Additional observations with ALMA in the coming years, providing a longer time baseline, are needed to confirm such a movement.

A circumplanetary disk around PDS 70 c. Isella et al. (2019) reported the detection of c_smm using ∼67 mas resolution observations. We confirm this detection with higher angular resolution observations that enable us to separate the emission from the outer disk. Given that the location of c_smm is very close to the existing Hα and near-IR measurements of PDS 70 c (Isella et al. 2019), and to the expected positions of PDS 70 c at the time of our observations (Figure 3), we interpret it as tracing the millimeter emission of dust grains located in a CPD. Assuming that c_smm is optically thin, its flux density can be converted into a dust mass estimate, for a given dust opacity and temperature. We note that if the emission is optically thick, such an assumption would provide a lower limit in the dust mass. The CPD temperature is also uncertain. It is determined by the sum of various sources of heating, namely viscous heating due to accretion of material through the CPD, accretion shocks, and external irradiation from both the planet and the star (Isella et al. 2014, 2019; Andrews 2021). Using 2 M_Jup, 2 R_Jup, and 1055 K as the mass, radius, and temperature of PDS 70 c (Wang et al. 2021), a mass accretion rate of 10⁻⁸ M_Jup/year (Haffert et al. 2019), we find that at a radial distance of 1 au from the planet, T_vis = 3 K, and T_p,irr = 18 K. Considering a stellar-irradiation temperature of T_s,irr = 24 K at the location of PDS 70 c (obtained from the radiative transfer model of Keppler et al. 2019), the CPD temperature at 1 au is ${T}_{\mathrm{CPD}}^{4}$ = ${T}_{\mathrm{vis}}^{4}$ + ${T}_{{\rm{p}},\mathrm{irr}}^{4}$ + ${T}_{{\rm{s}},\mathrm{irr}}^{4}$, that is T_CPD ∼ 26 K. Considering a typical dust opacity for 1 mm sized grains of 3.63 cm² g⁻¹ (Birnstiel et al. 2018) and a temperature of 26 K, we estimate a CPD dust mass of ∼0.007 M_⊕. A lower dust mass would be inferred if the dust temperature is higher than considered here (Schulik et al. 2020).

However, PDS 70 c is massive enough to carve a gap, and, as a consequence, large grains are trapped in a pressure maximum in the outer disk while small grains, well coupled to the gas, can flow inward. This is confirmed by the different cavity outer radii measured in scattered light compared to mm wavelengths (probing small and large grains, respectively; Keppler et al. 2019). Therefore, the CPD is only replenished with small dust particles that leak into the cavity (Bae et al. 2019) through meridional flows from the upper protoplanetary disk layers (e.g., Kley et al. 2001; Ayliffe & Bate 2009). If the CPD contains only small 1 μm sized grains (with an opacity of 0.79 cm² g⁻¹; Birnstiel et al. 2018) the CPD dust mass increases to ∼0.031 M_⊕. It is of course possible that the CPD hosts a range of particle sizes if the grains can grow. Bae et al. (2019) found that, if a steady state is achieved between the mass inflow to the CPD and the mass accretion rate onto the planet, the amount of sub-micron grains in the CPD would largely underestimate the observed mm flux and that accumulation of grains beyond the steady-state amount and/or in situ grain growth is needed to account for it. In Appendix D, we show the range of dust masses that the CPD would have for various dust grain size distributions, as a function of the maximum grain size. With these mass estimates, the ratio between the CPD dust mass and the planet mass, considering 2 M_Jup (Wang et al. 2021), ranges between 1 and 5 × 10⁻⁵.

If small grains can grow to mm sizes within the CPD, they could rapidly be lost as they efficiently drift toward the planet, and it only takes 100–1000 yr for an accreting CPD to lose all its mm dust (Zhu et al. 2018). However, as in protoplanetary disks, local gas pressure maxima can act as particle traps and prevent these grains from drifting. Interestingly, this can occur naturally in CPDs. Most of the gas that is feeding the CPD through meridional flows is then radially flowing outward in a decretion disk. The balance between the sub-Keplerian headwind and viscous outflow associated with a decretion flow leads to a global dust trap (Batygin & Morbidelli 2020). As a consequence, dust grains with sizes 0.1–10 mm may be trapped in the CPD and as the dust-to-gas ratio increases, streaming instabilities might be triggered (Drążkowska & Szulágyi 2018), or gravitational fragmentation in the outer regions of the CPD (Batygin & Morbidelli 2020) that will eventually lead to the formation of satellitesimals. At the same time, dust particles can accrete via pebble accretion onto the satellitesimals formed in situ or captured from the disk edge (e.g., Ronnet & Johansen 2020).

Our observations also put a strong constraint on the spatial extent of the CPD as seen in the dust emission at mm wavelengths. The emission c_smm is unresolved even at our highest angular resolution, and its peak intensity is similar over a range of beam sizes, until ∼40 mas, beyond which the CPD does not separate from the outer disk anymore (Appendix A.3). This indicates that it is more compact than 1.2 au in radius. On the other hand, there is a lower limit to the CPD extent needed to account for the observed flux. Assuming that it is a uniform disk with an optical depth of 1, and considering a temperature of 26 K, we find that it has a radius of 0.58 au. These two values (0.58 and 1.2 au) are therefore the lower and upper limits on the CPD radial extent constrained from our observations. The CPD is expected to be truncated (in gas) at one-third of the Hill radius, which for PDS 70 c, assuming a planet mass of 2 M_Jup at 34 au, is 1/3 × 3.1 ∼1 au. Three-dimensional (3D) simulations show that isothermal CPD are bound within 10% of the Bondi radius (Fung et al. 2019), that is 1/10 × 11 ∼1.1 au for PDS 70 c assuming a local temperature of 26 K. Both estimates, therefore, are consistent with our constraints. However, we cannot rule out the possibility that the gas component of the CPD extends beyond the dust component, in particular if some dust grains in the CPD drift inward.

Extended faint emission near PDS 70 b. The nature of the material close to PDS 70 b is unclear. It is not detected in the images obtained at high resolution with small synthesized beams, but is apparent at low S/N at intermediate resolution, indicating that it has a low surface brightness. It is confirmed in the two epochs 2017 and 2019, when combined with the short baselines data. b_smm appears offset toward the southwest from the position of PDS 70 b, confirming the findings of Isella et al. (2019), who speculated that it could be tracing dust trapped at the Lagrangian point L5 (Montesinos et al. 2020), if the planet is on an inclined orbit. The shape of the b_smm in our images suggests that it could also trace the faint signature of a streamer connecting the planets to the inner disk. Evidence for dust grains in the vicinity of PDS 70 b is clear already from the IR spectral energy distribution (Stolker et al. 2020; Wang et al. 2021), likely explaining the non-detection of Brγ (Christiaens et al. 2019) and Hβ emission lines (Hashimoto et al. 2020). It is interesting to understand why PDS 70 b, at the sensitivity of our observations, does not seem to host a compact, dusty, circumplanetary disk as PDS 70 c. A possibility would be that PDS 70 b has a much smaller Hill radius than PDS 70 c, as it orbits at smaller separation. Another natural explanation could be that PDS 70 b is starved of dust grains, as only the small grains that leak through the orbit of PDS 70 c and are transported through a streamer from the outer to the inner planet would enter the region of influence of PDS 70 b. Finally, it could be that the nature of the CPD is different around the two planets, with a decretion disk around PDS 70 c, and an accretion disk around PDS 70 b that is fed through a streamer coming from PDS 70 c rather than through meridional flows. More theoretical work looking at formation of CPDs in systems hosting two giant planets is needed to assess the potential differences between CPD formation in the inner and outer planet.

Inner disk. An inner dusty disk, evidenced in the IR spectral energy distribution and scattered light images is also clearly detected in our images up to ∼016 (∼18 au; see also Long et al. 2018; Keppler et al. 2019). Considering that the planets are filtering material from the outer disk such that only small dust particles can flow in the cavity, as for the CPD, it is unclear whether the inner disk mm emission is due to a population of small or large dust grains. To address this question, we computed the dust surface density and optical depth radial profiles of the continuum emission, using the combined data set (SB16+IB17+LB19) imaged with robust = 1. We consider four models for the dust grain population, that follow a size distribution n(a)da ∝ a^−3.5 da with a maximum grain size ${a}_{\max }$ of 10 μm, 100 μm, 1 mm and 1 cm, and a minimum size of 0.05 μm. We use the DSHARP opacities (Birnstiel et al. 2018) and the temperature profile output of the radiative transfer model of Keppler et al. (2018). The dust surface density as well as the total optical depth τ_ν is numerically computed, considering scattering and absorption opacities (Sierra & Lizano 2020; Sierra et al. 2021). The left panel of Figure 4 shows the total optical depth τ_ν for all four models. The right panel shows the dust surface density profiles (top) and corresponding cumulative masses (bottom). The dust surface density is maximum at the outer disk that is obviously the disk region that contributes to most of the dust mass (∼0.24 × 10⁻³ M_⊙ for ${a}_{\max }=1$ mm). We note that without the inclusion of scattering, the optical depth would follow the curve of the dust population with ${a}_{\max }=10$ μm, as the albedo at mm wavelengths is negligible for these small grains. In all these models, the inner disk is optically thin, with a total dust mass of ∼2 × 10⁻⁷–10⁻⁶ M_⊙ (i.e., 0.08–0.36M_⊕).

Zoom In Zoom Out Reset image size

Therefore, it appears that the emission at 855 μm from the inner disk regions located within the orbit of PDS70 b could be accounted for by a population of small grains. Interestingly, we note that the near-IR excess apparent in the spectral energy distribution of PDS70 is very low (Dong et al. 2012). This emission is mostly due to the thermal emission of small grains located within the innermost au, and such a low excess could indicate a low small-dust mass content in the inner disk and therefore suggest the additional presence of larger dust grains in order to account for the measured flux at 855 μm. However, the inner disk emission in the IR could still be optically thick (Dong et al. 2012), making it difficult to directly relate to our submillimeter observations and multiple wavelength observations in the millimeter regime are needed to constrain the grain size population in the inner disk. We note that the brightness temperature might be underestimated near the star because of our limited angular resolution and that it is possible that the innermost disk regions are optically thick also at submillimeter wavelengths. The longevity of the inner disk remains unclear; the replenishment flow is controlled by the planets, and if it is so strongly depleted (in gas) it may not allow grains to grow efficiently. It is possible that some of the dust in the inner disk is of a second generation produced by collision of larger bodies, perhaps stirred up by PDS 70 b. The star exhibits a small, but non-negligible, mass accretion rate, for which an additional mass reservoir in the inner disk, such as a dead zone, was recently suggested (Thanathibodee et al. 2020). Determining the physical conditions therein, in particular the dust-to-gas ratio, would be crucial for understanding whether such an inner disk can still grow terrestrial planets within a system hosting two outer giant planets. The current dust mass estimates are so low that it is unlikely that planets could form through pebble accretion (Lambrechts et al. 2019).

Outer disk structure. Our observations at high angular resolution indicate that the outer disk hosts substructures. In addition to an "arc" in the northwest, already seen at lower-resolution images (Long et al. 2018; Keppler et al. 2019), it resolves into two components, that can be either a double-ring structure with a dip at ∼055 or a bright ring with an inner shoulder. Interestingly, Huang et al. (2020) also found with high-resolution observations, a two-component structure in GM Aur, with a bright ring and an outer shoulder. It is unclear if such two-component structure in PDS 70 could be due to a secondary gap induced by PDS 70 c as an outer secondary gap opens only when the disk is sufficiently cold (Bae & Zhu 2018), with (h/r)_p ≲ 0.06 where (h/r)_p is the disk aspect ratio at the location of the planet ((h/r)_p ≃ 0.08 at PDS 70 c's location; Bae et al. 2019). On the other hand, recent 3D planet–disk interaction simulations including both gas and dust components showed that dust grains at the gap edge can have radial structures (Bi et al. 2021), potentially induced by corrugated vertical flows driven by the spiral wave instability (Bae et al. 2016a, 2016b) or meridional flows (Fung & Chiang 2016). Alternatively, such a substructure could be due to the presence of an additional, yet-undetected low-mass planet embedded within the outer disk. Similar multiple-ring substructures were also observed in other transition disks, such as HD 169142, in which three narrow rings were found and interpreted as tracing a migrating 10 M_⊕ in a low-viscosity disk (Pérez et al. 2019). However, hydrodynamical simulations show that thermodynamics can dramatically affect the structure of gas and dust, with different disk-cooling timescales leading to different planet-induced substructures (Facchini et al. 2020). Further chemical surveys will help to constrain the density and temperature structures (Facchini et al. 2021), enabling to test the possibility that an additional, low-mass planet is responsible for the structured outer ring and constrain the mass and radial location of that planet. We note that it is unlikely that an additional planet within the outer continuum ring disrupt the planetary system. In a two-planet system neglecting the eccentricity damping from the protoplanetary disk gas, the planets can avoid close encounters and are Hill-stable when their orbital separation is greater than 3.46 R_H, where ${R}_{{\rm{H}}}={a}_{1}{[({M}_{1}+{M}_{2})/3{M}_{* }]}^{1/3}$ is the mutual Hill radius (Gladman 1993; Barnes & Greenberg 2006). The addition of a third planet generally makes the stability criteria more stringent because the conservation of the total angular momentum and energy can no longer guarantee the avoidance of close encounters even for initially large separations beyond the Hill-stability criteria (Tamayo et al. 2015). However, provided that the protoplanetary disk gas provides sufficient eccentricity damping, Tamayo et al. (2015) argued that the two-planet criteria can still be used in three-planet cases. Assuming a range of 1–10 M_Jup for PDS 70 c and a Saturn mass for the hypothesized additional planet, this criteria is met when the latter is located beyond 44–53 au. Therefore, the system would be dynamically stable if the additional planet is located within the dip in the outer continuum ring at ∼60 au. Future numerical simulations will allow to further test our conclusions.

In this Letter, we report new ALMA observations obtained at high angular resolution (∼20 mas) at 855 μm of the PDS 70 system. We confirm the tentative detection by Isella et al. (2019) of a compact source co-located with the position of PDS 70 c with an independent data set at higher angular resolution. These new observations provide the most compelling evidence of the presence of a CPD around an accreting planet to date. Future molecular line IR observations at very high angular resolution may be able to detect rotating gas around PDS 70 c, providing conclusive results on the nature of the continuum mm emission. The detection of unresolved (r < 1.2 au) emission around planet c confirms that circumplanetary material is able to retain dust for long timescales, as required in satellite formation models.

These ALMA observations shed new light on the origin of the mm emission close to planet b. The emission is diffuse with a low surface brightness and is suggestive of a streamer of material connecting the planets to the inner disk, providing insights into the transport of material through a cavity generated by two massive planets. The non-detection of a point source around PDS 70 b indicates a smaller and/or less-massive CPD around planet b as compared to planet c, due to the filtering of dust grains by planet c preventing large amounts of dust to leak through the cavity, or that the nature of the two CPDs differ. We also detect a faint inner disk emission that could be reproduced with small 1 μm dust grains, and resolve the outer disk into two substructures (a bright ring and an inner shoulder).

PDS 70 is the best system to date to study and characterize circumplanetary disks, but also planet–disk interactions and disk cavity clearing by massive planets. The two massive planets, likely migrating outward in a grand tack-like scenario (Bae et al. 2019), are reminiscent of the Jupiter–Saturn pair, at larger distances from the star. Detailed studies of the circumplanetary disks, and of the leakage of material through the cavity, will provide strong constraints on the formation of satellites around gas giants, and on the ability to provide the mass reservoir needed to form terrestrial planets in the inner regions of the disk. Upcoming studies of the gas kinematics and chemistry of PDS 70 will complement the view provided by this work, serving as a benchmark for models of satellite formation, planet–disk interactions, and delivery of chemically enriched material to planetary atmospheres.

We are very grateful to Marco Tazzari and Antonella Natta for their support during the Covid-19 pandemic. We thank the North American ALMA ARC for their help. We acknowledge Ryan Loomis, Jason Wang, and Rens Waters for insightful discussions and the reviewer for helpful suggestions that improved the manuscript. This Letter makes use of the following ALMA data: ADS/JAO.ALMA#2018.A.00030.S., ADS/JAO.ALMA#2017.A.00006.S, ADS/JAO.ALMA#2015.1.00888.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 101002188 and No. 832428). J.B. acknowledges support by NASA through the NASA Hubble Fellowship grant #HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. S.F. acknowledges an ESO Fellowship. S.A. acknowledges support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. A.I. acknowledges support from the National Science Foundation under grant No. AST-1715719 and from NASA under grant No. 80NSSC18K0828. J.M.C. acknowledges support from the National Aeronautics and Space Administration under grant No. 15XRP15_20140 issued through the Exoplanets Research Program. N.T.K. and P.P. acknowledge support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research.

Facilities: ALMA - Atacama Large Millimeter Array, VLT(SPHERE) - , VLT(MUSE). -

Software: GALARIO (Tazzari et al. 2018), CASA (McMullin et al. 2007), Matplotlib (Hunter 2007), numpy (van der Walt et al. 2011), emcee (Foreman-Mackey et al. 2013), frank (Jennings et al. 2020).

A.1. Observing Log of the New ALMA Observations

A.2. Image Galleries and Corresponding Fluxes

To test the effect of the angular resolution on the image features and assess whether they are recovered in various images, we performed a grid of CLEANed and JvM-corrected images using different data sets, and different Briggs robust parameters. Figure 5 presents the resulting images. Corresponding image properties and fluxes are reported in Tables 4 and 5. Figure 6 shows the residual images (called "cavity images") obtained after subtracting the Fourier transform of the CLEAN model of the outer disk, for robust values of 0.5, 1, and 2.0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Summary of Disk and CPD Properties for Various Data Sets

Obs.				Disk		Emission around PDS70 c
Data Set	Briggs	Beam, Position Angle	rms Noise	Peak Iν	Total Flux	Peak Iν	Gauss Fit Peak Iν
	par.	(mas × mas)	(μJy beam−1)	(m Jy beam−1 )	(mJy)	(μJy beam−1)	(μJy beam−1)
LB19	0	22 × 22, 29°	20.4	0.29	172 ± 17	90 ± 20	91 ± 10
	0.3	26 × 25, 31°	14.6	0.33	175 ± 17	81 ± 15	82 ± 6
	0.5	29 × 27, 41°	11.0	0.37	176 ± 18	71 ± 11	71 ± 4
	1	42 × 34, 47°	8.2	0.65	196 ± 19	57 ± 8	49 ± 12
	2	47 × 40, 63°	6.2	0.82	193 ± 19	71 ± 6	60 ± 8
LB19+SB16	−0.5	20 × 20, 26°	26.3	0.27	184 ± 18	95 ± 26	97 ± 12
	−0.3	21 × 21, 2°	22.1	0.27	186 ± 19	88 ± 22	90 ± 11
	0	24 × 23, 30°	15.7	0.30	188 ± 19	86 ± 16	89 ± 8
	0.3	29 × 26, 40°	10.1	0.37	189 ± 19	82 ± 10	84 ± 6
	0.5	36 × 30, 44°	8.8	0.49	190 ± 19	80 ± 9	80 ± 8
	1	51 × 44, 63°	4.8	0.96	189 ± 19	81 ± 5	/
	2	60 × 54, 96°	3.9	1.37	189 ± 19	189 ± 4	/
IB17+SB16	−0.5	59 × 43, 59°	24.7	1.01	176 ± 18	105 ± 25	111 ± 25
	−0.3	60 × 44, 59°	20.5	1.06	176 ± 18	91 ± 20	98 ± 23
	0	64 × 48, 61°	15.5	1.20	176 ± 18	107 ± 15	100 ± 28
	0.3	70 × 54, 63°	11.0	1.45	176 ± 18	178 ± 11	182 ± 34
	0.5	75 × 59, 64°	9.1	1.68	177 ± 18	264 ± 9	428 ± 38
	1	87 × 69, 66°	6.3	2.22	178 ± 18	519 ± 6	714 ± 43
	2	93 × 74, 67°	5.0	2.48	178 ± 18	683 ± 5	817 ± 50
LB19+IB17+SB16	−0.5	24 × 23, 45°	16.7	0.29	173 ± 17	95 ± 17	95 ± 8
	−0.3	26 × 24, 41°	12.9	0.30	173 ± 17	89 ± 13	87 ± 6
	0	31 × 26, 44°	10.1	0.36	173 ± 17	86 ± 10	79 ± 8
	0.3	40 × 32, 47°	8.4	0.53	174 ± 17	83 ± 8	73 ± 9
	0.5	45 × 37, 51°	6.6	0.68	174 ± 17	79 ± 7	67 ± 9
	1	63 × 54, 78°	4.4	1.33	176 ± 18	170 ± 4	/
	2	70 × 63, 81°	3.5	1.68	176 ± 18	257 ± 3	/

Note. The c_smm peak intensities were computed with the CASA task imstat in an aperture centered on the CPD, with major/minor axis twice the beam major/minor axis (column 7) and with a Gaussian fit when possible (column 8). The rms is computed considering an annulus between 24 and 6''. We considered 10% calibration uncertainty as the flux uncertainty.

Download table as:  ASCIITypeset image

Table 5. Extended Flux in Cavity, from r = 1 Images

Obs.	Inner Disk			Emission around PDS70 b
Data Set	Major/Minor Axis, Position Angle	Peak Iν	Integrated Flux	Peak Iν	Integrated Flux
	(mas × mas)	(μJy beam−1 )	(μJy)	(μJy beam−1 )	(μJy)
LB19	129 ± 11/93 ± 8, 148 ± 11°	75 ± 6	719 ± 60	/	/
LB19+SB16	128 ± 11/94 ± 9, 152 ± 12°	126 ± 9	817 ± 69	46 ± 5	101 ± 10
IB17+SB16	102 ± 12/81 ± 13, 171 ± 32°	367 ± 18	888 ± 59	56 ± 6	40 ± 4
LB19+IB17+SB16	117 ± 5/91 ± 4, 166 ± 8°	174 ± 5	727 ± 27	27 ± 4	38 ± 3

Note. The inner disk properties were derived using a Gaussian fit in an elliptical mask centered in the central pixel, sized 015 × 012. Deconvolved major, minor axis FWHM, and position angle are given. For the material around PDS70 b, we defined a rectangular aperture of 008 × 015, with PA = 55°.

Download table as:  ASCIITypeset image

A.3. Peak Intensity of the Continuum Emission Associated with the Planets

Figure 7 shows the peak intensity of c_smm as a function of angular resolution. Depending on the data set, and the robust parameter, its peak intensity ranges between 80 ± 6 and 107 ± 15 μJy beam⁻¹ when it is well separated from the outer ring. At larger resolution than ∼60 mas, the peak intensity increases because the beam contains contribution from the outer disk. The gray area reports the estimate of Isella et al. (2019). In contrast the peak intensity of b_smm varies between 46 ± 5, 56 ± 6, and 27 ± 4 μJy beam⁻¹ for three different data sets (LB19+SB16, IB17+SB16, and LB19+IB17+SB16, respectively) imaged at resolutions of 51 mas × 44 mas, 87 mas × 69 mas, and 63 mas × 54 mas, respectively (see Table 5).

Zoom In Zoom Out Reset image size

To better assess the presence of substructures in the outer disk, we fit azimuthally averaged deprojected visibilities using the frank package that models an axisymmetric surface density profile using a flexible Gaussian process (Jennings et al. 2020). To do so, we considered the combined data set LB19+IB17+SB16, which has the best uv plane coverage, assuming that the outer disk brightness distribution has not changed between these observations. The data was averaged into intervals of 30 s and 1 channel per spectral window to reduce data volume.

As the disk presents an asymmetric arc-like feature in the northwest that can lead to overestimate the disk radial intensity profile when fitting with an axisymmetric model, we followed the procedure developed in Andrews (2021) that modifies the visibilities by removing a model for the asymmetry before fitting with frank. We mask the emission between position angles of −85° and 40°. The asymmetry is selected in the CLEAN model image, and the mean radial profile of the CLEAN model from the disk outside the asymmetric region is subtracted from the model image, allowing us to obtain a model for the asymmetry only. The Fourier transform of the asymmetry model was then subtracted from the original visibilities and the resulting set of visibilities are further modeled. frank fits deprojected visibilities, and inaccurate estimates of the geometric parameters for the deprojection, (ΔR.A., Δdecl., inc, PA), can lead to significant residuals. Automatic procedures performed poorly to find the best parameter set, and therefore we optimized those parameters by hand, exploring different values of spatial offset and geometry in steps of 1 mas and 05. The final values adopted for the frank fit are (12 mas, 15 mas) for (dRA, dDec), and (49.5, 161.0) for (inc, pa). We tested the sensitivity of the fit to the hyperparameters α (varied between 1.05 and 2.00) and w_smooth (varied between 10⁻⁴ to 10⁻¹) and found no significant difference on the residuals. Therefore, we fixed for standard values w_smooth = 10⁻⁴ and α = 1.30.

The fit to the data and the corresponding profile are shown in Figure 8 (top row). The best-fit model indicates a radial profile with two local maxima for the emission of the outer disk, confirming the findings of Keppler et al. (2019) with lower-resolution observations. At the angular resolution of our observations, the two peaks are separated by ∼7 beams. However, it is unclear whether the outer disk is composed of two separated broad rings, or of a ring with an inner shoulder. No clear gap or ring is evidenced in the inner disk. We note the presence of a possible additional shoulder at 085. The model and residuals are imaged exactly as the observations, with a robust parameter of 0.5, and are presented in the bottom row of Figure 8. The residuals indicate that the axisymmetric model does not account for the full complexity of the emission, in particular the disk can be vertically thick and flared. A dedicated two-dimensional modeling of the visibilities, which is beyond the scope of this Letter, is needed to properly assess the morphology of the disk and will be presented in a forthcoming study.

Zoom In Zoom Out Reset image size

Table 6 provides the published astrometry of PDS 70 b and PDS 70 c, as well as the predicted locations at the time of our observations based on the best orbital fits by Wang et al. (2021).

We consider three models for the dust grain population in the CPD around PDS 70 c, that follow different size distribution n(a)da ∝ a^−3.5 da, ∝ a^−3.0 da, and ∝ a^−2.5 da and show the predicted dust mass as a function of the maximum grain size ${a}_{\max }$ in Figure 9. We consider a minimum grain size of 0.05 μm, a temperature of 26 K and use the DSHARP opacities (Birnstiel et al. 2018).

Zoom In Zoom Out Reset image size